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Abstract 

A 1-d Ising model is shown to reproduce qualitatively the dynam- 
ics of ripple formation. Saltation effect is imposed using a Kawasaki 
dynamics and a pair interaction over some distance £. Within this 
model, the ripple state turns out to be metastable in agreement with 
cellular automata simulations as well as recent under water experi- 
ments. A dynamical phase diagram is obtained. A mean-field solution 
of the problem is given in terms of the ripple size. A mapping is then 
performed onto a two-dimensional i x oo static problem. 
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1 Introduction 



Ising spin models have been used to describe a tremendous spectrum of phys- 
ical and non-physical problems [p], [|. Among others, they have proven useful 
to understand some aspect of granular media. In particular compaction dy- 
namics under tapping |L]. 0. 

Herein, we present a simple 1-d Ising model to describe qualitative fea- 
tures of ripple formation. Previous experimental and theoretical works have 
underlined the primary role played by saltation in the emergence of sand rip- 
ples and dunes [§, ||, ||, |7| . Saltation means grain ejection and eolian transport 
over long distances. The hopping and rolling of grains have generally led to 
travelling ripple structures which in turn merge and grow. 

In the following, an Ising spin description of the dynamics of ripples 
formation is presented. A 1-d Hamiltonain is suggested. It is studied numer- 
ically with a Kawasaki dynamics. A mean-field treatment is also performed. 
A dynamical phase diagram is obtained. Results are discussed and com- 
pared with recent cellular automaton results M and with recent underwater 
experiments Q. 

2 The model 

A natural way to map a one-dimensional granular landscape into a spin model 
is to consider a 1-d Ising model {crj = ±1}. Let the wind blowing from left 
to right. At each site i, <ii = +1 if the local slope is positive (exposed to 
the wind). It is <7j = —1 if the local slope is negative (screened by the ripple 
crest and unaffected by the wind). Successive slopes being then associated 
to the succession of spin domains. Spin domains represent the ripple sides. 
A sketch of such a mapping is shown in Figure 1. 

In Figure 1, a grain is first shifted over a distance i = 7 by saltation. An 
avalanche then occurs to relaxe the surface. From the magnetic viewpoint, 
both extraction and relaxation of one grain produce two pairs of spin flips 
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meaning that the whole dynamics goes by spin pairs. Moreover, it proves that 
the dynamics does not modify the total magnetization m = J^i^i- To embody 
these physical effects charaterizing the saltation-avalanche mechanism, we 
consider the so-called Kawasaki dynamics || . Having a 1-d Ising system the 
equilibrium total magnetization is zero. In addition, it will stays zero since 
the Kawasaki dynamics preserves the total magnetization. 

Avalanches which relaxe the granular surface can be viewed as small ther- 
mal fluctuations along the ripple sides. In our spin mapping such a scheme 
corresponds to a ferromagnetic nearest neigbor coupling J nn divided by ksT 
where T is the temperature and ks the Boltzman's constant. Moreover, 
saltation occurs only on the faces exposed to the wind and leads to a per- 
turbation at the position i + i where I is the distance overwhich the sand 
is transported by saltation. As for avalanches, we mimic this effect by a 
ferromagnetic coupling J s which connects a site i to a site i + £, divided by 
k B T. 

On this basis we consider the following dimensionless Hamiltonian 

+00 +00 
E = -K nn - K s ^2 ai(j i+e 6 a . }+ i , (1) 

i=— oo i=— oo 

where 5 x>y is the Kronecker function, K nn = and K s = The Kro- 
necker delta function in the second term assumes that only faces exposed to 
the wind (spins +1 for a wind blowing from left to right) play a role for salta- 
tion. These assumptions seem reasonable and not of drastic consequences at 
our present level of investigation. 

3 Numerical simulations 

Numerical simulations have been performed on a 1-d system with periodic 
boundary conditions, i.e. a loop. The initial spin loop configuration is ran- 
dom with m — 0. Thereafter, one positive spin and one negative spin are 
randomly selected at each time step. The energy variation AE due to the 
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flipping of both selected spins is calculated. Periodic boundary conditions 
are used in the energy calculation. The selected spins are then flipped with 
some probability p = ^^fg as in classical Monte Carlo methods || . The 
fact that two opposite spins are chosen at each time step preserves the m = 
value of the order parameter during the whole simulation, i.e. a Kawasaki 
dynamics. 

Since the magnetization m is zero and constant, a dynamical "order pa- 
rameter" can be defined with the domain wall density p, i.e. the density of 
nearest neighboring antagonist pairs of spins. One should note that 

where £ is the mean size of the spin domains. Loop lengths up to L = 
10000 spins have been used. Our numerical results has been checked against 
variation in the parameters (K nn , K s and £). Finite-size effects as well as the 
spin pattern time evolution have also been analyzed, it is worth noting that 
I and L have been always chosen to be incommensurate in order to avoid 
any artefact due to the periodic boundary conditions. 

Figure 2 presents typical space-time patterns in which domain walls are 
drawn for characterisitic lengths L = 144 and I = 21. Three different sets 
of parameter values K nn and K s are shown. Depending on the parameter 
values three different regimes are obtained. Either disorder (Figure 2a), or 
order with domains of roughly the size £/2 (Figure 2b), or also order with 
larger domains (Figure 2c). Domains having a size of the order £/2 are 
called thereby "ripples" . In parallel much larger domains are called "dunes" 
in reference to simulations in ||. 

Figure 3 contains the domain wall density p as a function of time for 
the three regimes of Figure 2. It is observed: (i) an erratic signal for the 
disordered regime presenting huge fluctuations of p, (ii) an apparently stable 
density p ~ t/2 after a short transient regime and (iii) a low density regime 
(p ~ 0) decorated with small fluctuations after some transient regime. 
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Figure 4 is the associated resulting dynamical phase diagram as function 
of (K nn , K s ) obtained numerically for £ > 1. Three different dynamical 
regimes are thus found. A region of the disorder regime is represented in 
grey in the phase diagram. The boundary with the order regime is found to 
be stable whatever the values of £ > 1 and L. Curves p(K nn ) are shown in 
Figure 5 for different values of K s where one can clearly see a crossover around 
K nn m 1. At this point, the domain size is equivalent to £/2. Moreover, this 
boundary is stable in time. This result could be astonishing since finite- 
size effects are expected. One should recall that m = corresponds to the 
equilibrium state of a one-dimensional Ising model and thus fluctuations and 
the domain wall evolution take place in the equilibrium state. Those are 
mainly driven by the Ising interaction K nn . 

For the order regime (white region in Figure 4), two distinct regions 
coexist: ripples (for K s » K nn ) and dunes (for K nn » K s ). The dashed 
line in Figure 4 represents the boundary between these regimes. However, 
this boundary moves slightly when the system size L or the saltation length 
£ is changed. In fact, the ripple state is metastable. After very long times, 
ripples can suddenly merge into dunes (as seen by geologists). The larger 
K s , the more stable ripples are. 

Different initial correlations have been examined. Starting from a perfect 
antiferromagnetic situation, the system reaches the same regimes as described 
above. The dynamics is similar to what we have observed for the random 
inital situation. However, if one starts with a system containing only two 
large domains, i.e. a single giant dune, the system never produces ripples 
but remains unchanged in the order region of the dynamical phase diagram. 
In the disorder regime (in grey in Figure 4), the initial giant dune is however 
destroyed. 

Finally, one should note the case £ = 1 is quite particular. In that case, 
spin domains are produced even when K nn = for which a disorder regime is 
expected. As soon as £ > 1, the results summarized in Figure 4 are recovered. 
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4 Mean-field treatment 



Eventhough the Hamiltonian (1) is rather simple, it does not allow for an 
easy exact solution. For instance, the classical method of the transfer matrix 
|T3| ] leads to searching for 2 £+1 x 2 e+1 matrices. In order to solve this problem, 
we have considered as a first approximation that the mean domain size £ is 
the relevant parameter. Fluctuations in £ are neglected. 

Assuming an order with domains of size £, the first product OiOi+x of 
Eq.(l) gives due to the presence of domain walls which are energetically 
defavorable. The second product (Jidi+i of Eq.(l) gives a periodically oscil- 
lating function of period 2£. Thus, the energy per lattice site E/L can be 
written as 

— « -K nn ( ^— - ) - ^cos f ^— J , (3) 



L ""V £ 

which is drawn in Figure 6 as a function of £ and for arbitrary values of the 
coupling energies K nn and K s . Several minima are observed. The selection 
of various lengths £/2, £/3, £/&,... is thus expected from the above shape of 
E(£). The more stable local minima is for 1/2, i.e. it corresponds to ripples 
of size I. 

For £ > £/2, a large potential barrier separates the ripple state from the 
"giant dune" state. The asymmetry of this barrier implies the irreversibility 
of the "ripple — > dune" process. The ripple metastability is thus also obtained 
from a mean-field treatment in agreement with the above numerical results. 

Nevertheless, the mean-field solution does not predict any disordered dy- 
namical regime. Indeed, fluctuations of £ are neglected. Fluctuations destroy 
the order when K q is small. 
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5 Mapping onto a two-dimensional £xoo static 
problem 

Although ripple formation is a dynamical phenomenon, our study is gener- 
ated via Monte Carlo (equilibrium-like) simulations. The dynamical aspect 
survives only because the relaxation towards equilibrium is slow and ham- 
pered by metastability. 

Another remark concerns Eq.(l). Upon setting <5 CTi ,+i = |(1 — <7j) in 
Eq.(l), multiplying out the a's, and using the fact that m = J2i< J i one finds 

E = -K nn <Ji<7i+i - K s ^ OiOi+f, + ^K s m (4) 

i i 

in which the last term is an irrelevant constant that is set equal to zero later. 
Replace now each site index i by two integer coordinates j and k uniquely 
defined by the conditions 

i = j£ + k with k e {0,1, ...,£- 1} (5) 

and write <Tj = a^ k . In terms of these new variables the Hamiltonian Eq.(4) 
takes the form 

oo t— 1 oo I— 1 

E = ~ K nn Vj,k<7j,k+1 ~ K s Yl Yl a j,k<7j+l,k (6) 

j=— oo k=0 j=—oo k=0 

provided we interpret (j, k) as a site on an £ x oo cylindrical lattice with 
screw boundary conditions (CL/sbc) around the cylinder. Hence E of Eq.(6) 
describes Onsager's anisotropic nearest-neighbor Ising model on a CL/sbc, 
but with the additional constraint m = 0. 

This sheds an interesting light on the simulation results. The relevant 
facts are 

(a) The constraint m = and the sbc are irrelevant for the thermody- 
namics of this model, which is therefore the same as that of the Onsager 
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model without constraint and in magnetic field zero. 



(b) The I x oo lattice has no phase transition. But in the limit I — > oo 
this model does have the ferromagnetic Onsager phase transition. Of course 
for I sufficiently large simulations will show an order/disorder transition-like 
behavior. 



(c) The ground state Q , due to the condition m — 0, necessarily has an 
interface. If the j direction (the one parallel to the cylinder axis) is infinite, 
Qq is a cylinder which is positively magnetized for j > and negatively 
magnetized for j < 0. If the j direction is finite (let's say —J < j < J 
for some J), then Q is still a low-lying state of energy E(Q ) ~ —£K S , 
but another low-lying state appears viz. the state Q\ which has positive 
magnetization for E/2 < k < i—1 and negative magnetization for < k < E/2 
(£ supposed even for convenience). Its energy is E{Qi) ~ —2JK nn . Both 
states are local energy minima in configuration space (because in both the 
interface length is at a local minimum), and, with the Kawasaki dynamics 
of this paper, the higher one of the two states is metastable with respect to 
the lower one. It is evident that Q Q corresponds to the " dune" and Q\ to the 
"rippled state" . 



6 Connexion to recent underwater experiments 

The physics underlying both ripple and dune formation is far from being 
understood making any analytic result valuable. Indeed, it has been believed 
that dune dynamics are of a different nature than ripple ones. Moreover the 
question of ripple growing to eventually produce a dune is still open. 

Bagnold [[10] distinguished ripples and ridges. Periodicity of the former 



ones was time independent while periodicity of the latter ones was growing 
with time. The present work suggests in a simple way that both patterns 
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are two aspects of the same phenomenon since the ripple state seems to be 
metastable. 

Recent experiments of underwater sand @, |Tl| have exhibited the forma- 
tion of ripples. There, sudden merging of ripples do also occur in qualitative 
agreement with our results (see the upper part of the dynamical phase dia- 
gram of Figure 4). Moreover, these authors have found that a "transition" 
takes place between a flat landscape regime and ripple formation at a critical 
value of the Reynolds number Re. This feature is to be put in parallel to 
the transition that we have obtained (see Figure 4) between a disordered 
regime and a ripple dynamical regime. Indeed, the saltation length i could 
be directly connected to the physical parameter Re. 



7 Comparaison with a cellular automaton 

In addition, a recent cellular automaton model || did produce a dynamical 
ripple state which eventually stabilizes into a giant dune. This model consid- 
ers a granular landscape for which surface local fluctuations are dissipated 
through both saltation and avalanches [I2|. Saltation is restricted to grains 
exposed to the wind while avalanches occurs everywhere. 

The associated dynamics is rather complex. First, ripples appear. Sec- 
ondly, they coalesce very slowly into "giant dunes" . Main features of ripples 
and dunes formation are recovered within the present spin framework except 
for the disordered regime which is absent in ||. 

Eventhough the cellular automaton model as well as actual spin-like 
model are far from real desert dunes and ripples, they seems to reproduce 
some analogous features. 
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8 Conclusion 



In summary, we have proposed a 1-d Ising model for ripple formation. The 
spin model has been studied numerically. In addition, a mean-field solution 
for the ripple size was also obtained. Moreover, the ripple state has been 
found to be metastable. A result which is in full agreement with very recent 
experiments under water as well as a previous landscape simulations. 

Nature is obviously much more complex that a one-dimensional spin 
model. However our very simple description reproduces some basic features 
of dune formation as seen above. Moreover, the model could be extended 
to include, for instance, a distribution of saltation lengths or even, different 
grain types. This is left for future work. We do hope our results will stimu- 
late additional experimental work in particular to check our dynamical phase 
diagram. 
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Figure Captions 



Figure 1 — (top) Sketch of the rule. A grain at position i is unstable and 
is carried over by the wind to the position i + £ where an avalanche relaxes 
the surface. Spin signs indicate the slope of the landscape, (bottom) After 
the saltation- avalanche event, the magnetization remains unchanged. 

Figure 2 — Typical space-time patterns in which domain walls are drawn: 
L = 144 and I = 21. Three different sets of parameter values are shown: (a) 
K nn = 0.5 and K s = 0.5, (b) K nn = 2.0 and K s = 2.5, and (c) K nn = 2.5 
and K s = 0.5. 

Figure 3 — The density of domain walls p as a function of time. Three 
different sets of parameter values are shown, from top to bottom: K nn = 0.5 
and K s = 0.5, K nn = 2.0 and K s = 2.5, and K nn = 2.5 and K s = 0.5. 

Figure 4 — The dynamical phase diagram (K nn , K s ) of the model. The 
disorder phase is represented in grey in the diagram. The order phase is 
represented in white. Ripples and dunes are also distinguished. 

Figure 5 — The density of domain walls p as a function of the Ising interac- 
tion K nn . Three different values of K s are shown. The crossover separating 
disorder and order regimes is emphasized. 

Figure 6 — Typical shape of the energy per lattice site E/L as a function 
of the normalized domain size 
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